Spatio-temporal dynamics of coupled array of Murali-Lakshmanan-Chua circuits 
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The circuit recently proposed by Murali, Lakshmanan and Chua (MLC) is one of the simplest non- 
autonomous nonlinear electronic circuits which shows a variety of dynamical phenomena including 
various bifurcations, chaos and so on. In this paper we study the spatio-temporal dynamics in one 
and two dimensional arrays of coupled MLC circuits both in the absence as well as in the presence of 
external periodic force. In the absence of any external force, the propagation phenomena of travelling 
wave front and its failure have been observed from numerical simulations. We have shown that the 
propagation of travelling wave front is due to the loss of stability of the steady states (stationary 
front) via subcritical bifurcation coupled with the presence of necessary basin of attraction of the 
steady states. We also study the effect of weak coupling on the propagation phenomenon in one 
dimensional array which results in the blocking of wave front due to the existence of a stationary 
front. Further we have observed the spontaneous formation of hexagonal patterns (with penta-hepta 
defects) due to Turing instability in the two dimensional array. We show that a transition from 
hexagonal to rhombic structures occur by the influence of external periodic force. We also show the 
transition from hexagons to rolls and hexagons to breathing (space time periodic oscillations) motion 
in the presence of external periodic force. We further analyse the spatio-temporal chaotic dynamics 
of the coupled MLC circuits (in one dimension) under the influence of external periodic forcing. 
Here we note that the dynamics is critically dependent on the system size. Above a threshold size, 
a suppression of spatio-temporal chaos occurs, leading to a space-time regular (periodic) pattern 
eventhough the single MLC circuit itself shows a chaotic behaviour. Below this critical size, however, 
a synchronization of spatiotemporal chaos is observed. 



I. INTRODUCTION 



Coupled nonlinear oscillators in the form of arrays arise 
in many branches of science. These coupled arrays are 
often used to model very many spatially extended dy- 
namical systems explaining reaction-diffusion processes. 
One can visualize such an array as an assembly of a num- 
ber of subsystems coupled to their neighbours, or as cou- 
pled nonlinear networks (CNNs) [Chua & Yang, 1988; 
Chua, 1995]. Theoretically one can consider a system 
of coupled ordinary differential equations (odes) to rep- 
resent a macroscopic system, in which each of the odes 
corresponds to the evolution of the individual subsys- 
tem. Specific examples are the coupled array of an- 
harmonic oscillators [Marin & Aubry, 1996], Josephson 
junctions [Watanabe et al., 1995], continuously stirred 
tank reactors (CSTR) exhibiting travelling waves [Dolnik 
& Marek, 1988; Laplante & Erneux, 1992], propagation 
of nerve impulses (action potential) along the neuronal 
axon[Scott, 1975], the propagation of cardiac action po- 
tential in the cardiac tissues [Allesie et al., 1978] and so 
on. For the past few years several investigations have 
been carried out to understand the spatio-temporal be- 
haviours of these coupled nonlinear oscillators and sys- 
tems. The studies on these systems include the travelling 
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wave phenomena, Turing patterns, spatio-temporal chaos 
and synchronization [Feingold et al, 1988; Munuzuri et 
al, 1993; Perez-Munuzuri et al., 1993; Sobrino et al, 
1993; Perez-Munuzuri et al., 1994, Perez-Munuzuri et 
al, 1995; Kocarev & Parlitz, 1996]. 

One of the interesting behaviours of these coupled os- 
cillators is that they can exhibit a failure mechanism for 
travelling waves which cannot be observed in the homoge- 
neous continuous models [Keener, 1987]. Similar failure 
phenomena in the nerve impulse propagation have been 
observed in biological experiments [Balke et al., 1988; 
Cole et at, 1988]. 

Recently various spatio-temporal patterns have been 
found in discrete realizations of the reaction-diffusion 
models by means of arrays of coupled nonlinear elec- 
tronic circuits [Perez-Munzuri et al., 1993]. For exam- 
ple Perez-Muhuzuri and coworkers studied the various 
spatio-temporal patterns in a model of discretely coupled 
Chua's circuits. This system is a three variable model and 
coupled to its neighbours by means of a linear resistor. In 
this paper we consider similar arrays of coupled nonlin- 
ear electronic circuits consisting of Murali-Lakshmanan- 
Chua (MLC) circuits as a basic unit, which is basically 
a two variable model. The dynamics of the driven MLC 
circuit including bifurcations, chaos, controlling of chaos 
and synchronization has been studied by Murali et al. 
[1994a; 1994b; 1995; Lakshmanan & Murali, 1996]. Of 
particular interest among coupled arrays is the study 
of diffusively coupled driven systems as they represent 
diverse topics like Faraday instability [Lioubashevski et 
al., 1996], granular hydrodynamics [Umbanhawar et al., 
1996; Kudrolli et al, 1997], self organized criticality [Bak 
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& Chen, 1991] and so on. Identification of localized struc- 
tures in these phenomena has been receiving considerable 
attention very recently [Fineberg, 1996]. Under these cir- 
cumstances, study of models such as the one considered 
in this paper can throw much light on the underlying 
phenomenon in diffusively coupled driven systems. 

In the following by considering one and two dimen- 
sional arrays of coupled MLC circuits, we have studied 
the various spatio-temporal patterns in the presence as 
well as in the absence of external force. We present a 
brief description on the one and two dimensional arrays 
of MLC circuits in Sec. HU Then, in Sec. Hill in the 
absence of periodic external force, we study the propaga- 
tion phenomenon of travelling wave front and its failure 
below a critical diffusive coupling coefficient. We show 
that the propagation phenomenon occurs due to a loss of 
stability of the steady states via subcritical bifurcation 
coupled with the presence of necessary basin of attrac- 
tion of the steady states for appropriate diffusive coupling 
coefficient. We also study the effect of weak coupling on 
the propagation phenomenon which causes a blocking in 
the propagation due to the presence of a stationary front 
at the weakly coupled cell. In addition we also study the 
onset of Turing instability leading to the spontaneous 
formation of hexagonal patterns in the absence of any 
external force. 

Further as mentioned earlier it is of great physical in- 
terest to study the dynamics of the coupled oscillators 
when individual oscillators are driven by external forces. 
In Sec. IIVI we study the spatio-temporal patterns in 
the presence of periodic external force and investigate 
the effect of it on the propagation phenomenon and Tur- 
ing patterns. Depending upon the choice of control pa- 
rameters, a transition from hexagons to regular rhombic 
structures, hexagons to rolls and then to breathing os- 
cillations from hexagons are observed. The presence of 
external force with sufficient strength removes the penta- 
hepta defect pair originally present in the spontaneously 
formed hexagonal patterns leading to the formation of 
regular rhombic structures. The inclusion of external pe- 
riodic force can also induce a transition from hexagons to 
rolls provided there are domains of small roll structures 
in the absence of force. We further show that in the re- 
gion of Hopf- Turing instability, the inclusion of external 
periodic force with sufficiently small amplitude induces a 
type of breathing oscillations though the system shows a 
regular hexagonal pattern in the absence of any external 
force. 

Finally, we also study the spatio-temporal chaotic dy- 
namics of the one dimensional array of MLC circuits in 
Sec. Ivl when individual oscillators oscillate chaotically. 
In this case, the emergence of spatio-temporal patterns 
depends on the system size. For larger size, above a crit- 
ical number of cells, we observe a controlled space-time 
regular pattern eventhough the single MLC circuit itself 
oscillates chaotically. However, synchronization occurs 
for a smaller system size, below the threshold limit. Fi- 
nally, in Sec. I VII we give our conclusions. 



II. ARRAYS OF 
MURALI-LAKSHMANAN-CHUA CIRCUITS 

As mentioned in the introduction the circuit proposed 
by Murali, Lakshmanan and Chua (MLC) is one of the 
simplest second order dissipative nonautonomous circuit 
having a single nonlinear clement [Murali et ai, 1994a]. 
Brief details of its dynamics are given in the Appendix 
A for convenience. Here we will consider one and two 
dimensional arrays of such MLC circuits, where the in- 
tercell couplings are effected by linear resistors. 



A. One-dimensional array 

Fig. Ola) shows a schematic representation of an one 
dimensional chain of resistively coupled MLC circuits. 
Extending the analysis of single MLC circuit as given 
in the Appendix A, the dynamics of the one dimensional 
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FIG. 1: (a) Circuit realization of an one dimensional ar- 
ray of coupled MLC circuits, (b) Graph of the function 
h(x) Vs. x for the parametric choice {mo,mi,m2,e} = 
{-1.02, -0.55, -0.55,0} in Eq. (2). (c) Graph showing h(x) 
Vs. x for the parametric choice {mo, mi, m,2, e} = {—2.25,1.5, 
0.25,0} in Eq. (2). 

chain can be easily shown to be governed by the following 
system of equations, in terms of the rescaled variables 
(see Appendix A), 

ii = yi-h(x i )+D(x i+1 + x i - 1 -2xi), (la) 
ili = -<7Vi - 0Xi + Fsinut, i= 1,2, ••• ,N, (lb) 

where D is the diffusion coefficient, N is the chain length 
and h(x) is a three segment piecewise linear function rep- 
resenting the current voltage characteristic of the Chua's 
diode, 

!e + TR2X + (mo — mi), x > x 2 
e + mox, X\ < x < X2 (2) 

e + m\x — (mo — mi), x <x\. 

In |J2J m , mi and m 2 are the three slopes. Depending 
on the choice of mo, mi and mi one can fix the charac- 
teristic curve of the Chua's diode. Here e corresponds to 
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the dc offset. Of special importance for the present study 
are the two forms of h(x) shown in Figs, ^b) and^c). 
Fig. IHb) corresponds to the choice of the parameters 
{mo, mi, ma, e} = {-1.02,-0.55,-0.55,0} while Fig. 
^c) corresponds to the parameters {mo,mi,m2,e} = 
{-2.25,1.5,0.25,0}. As noted in the Appendix A, the 
significance of the characteristic curve in Fig. ^b) is 
that it admits interesting bifurcations and chaos, while 
Fig. corresponds to the possibility of bistable states 
with asymmetric nature which is a prerequisite for ob- 
serving wave front propagation and its failure. In the 
following studies we will consider a chain of N = 100 



MLC circuits. 



B. Two-dimensional array 

As in the case of one dimensional array one can also 
consider a two dimensional array with each cell in the 
array being coupled to four of its nearest neighbours with 
linear resistors. The model equation can be now written 
in dimcnsionlcss form as 



-h3 



yij - h(xij) + Di(x i+ ij + Xi-ij + x hj+ i + . 



'hj—i 



= -o"!Ji,j ~ Pxi.j + D 2 (y t +ij + Vi-i,j + Vi,j+i + Ut.j-i - + Fsxn.u)t = g(x it j,y it j), 

I 



i,i = 1,2, • 



(3a) 
,7V. 



This two dimensional array has N x N cells arranged in 
a square lattice. In our numerical study we will again 
take N = 100. 

In the following sections we present some of the inter- 
esting dynamical features of the array of MLC circuits 
such as propagation of wave front and its failure, effect of 
weak coupling in the propagation, Turing patterns, effect 
of external periodic forcing in the Turing patterns and 
spatio-temporal chaotic dynamics. We have used zero 
flux boundary conditions for the study of propagation 
phenomenon and Turing patterns and periodic boundary 
conditions for the study of spatio-temporal chaos. 



III. SPATIO-TEMPORAL PATTERNS IN THE 
ABSENCE OF EXTERNAL FORCE: 
TRAVELLING WAVE PHENOMENON AND 
TURING PATTERNS 

Transport processes in living tissues, chemical and 
physiological systems have been found to be associated 
with special types of waves called travelling waves. Ear- 
lier, continuous models were created to describe the wave 
propagation phenomenon in these systems, but these 
failed to cover all the important aspects of travelling 
wave behaviour. One of the most important of them 
is the so called travelling wave propagation failure, oc- 
curring at weak coupling between cells. It has been 
proved by Keener [1987] that propagation failure cannot 
be observed in a continuous, one variable, homogeneous 
reaction-diffusion system. Therefore to study this kind of 
phenomenon, one has to use discrete models. Recently 
travelling wave phenomenon and its failure have been 
studied in arrays of discretely coupled Chua's circuits 
[Perez-Muhuzuri et at, 1993]. Wave propagation and its 
failure have also been observed even in one variable mod- 
els like discrete Nagumo equation [Leenaerts, 1997]. In 



this section we carry out a study of such wave propa- 
gation phenomenon in one and two dimensional arrays 
of Murali-Lakshmanan-Chua's (MLC) circuit oscillators, 
without forcing, and investigate the mechanism by which 
such a failure occurs. 



A. Propagation phenomenon and its failure in 
one-dimensional array 

In order to observe wave front propagation in the one 
dimensional array of MLC circuits, we numerically in- 
tegrated the Eqs. (1) using fourth order Runge-Kutta 
method. In this analysis we fix the parameters at 
{j3,a,m ,m 1 ,m 2 ,e, F} = {1.0, 1.0, -2.25, 1.5, .25, 0, 0} 
so that the system admits bistability and this choice leads 
to the asymmetry characteristic curve for the Chua's 
diode as shown in Fig. ^c). The existence of bistable 
state in the asymmetric case is necessary to observe a 
wave front. Zero flux boundary conditions were used in 
the numerical computations, which in this context mean 
setting xq = X\ and xn + i = x^ at each integration step; 
similar choice has been made for the variable y also. To 
start with, we will study in this section the force-free 
case, F = 0, and extend our studies to the forcing case 
(F ^ 0) in the next section. 

The choice of the values of the parameters guaran- 
tees the existence of two stable equilibrium points = 
{ct(toi - mo — e)/(/3 + 77120"), /3(to -toi -e)/(/? + TO 2 0)} 
and Pf — {o"(toi — mo — e)/(j3 + mio"), /3(mo — mi — 
e)/(/3 + TOiO")} for each cell. In the particular case, corre- 
sponding to the above parametric choice, each cell in the 
array has three equilibrium points at = (3.0,-3.0), 
Pr = (-1.5,1.5) and if = (0,0). Out of these three 
equilibrium points, P^~ and P~ are stable while P® is 
unstable. Due to the asymmetry in the function h(x) for 
the present parametric choices, the basin of attraction of 
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the point P^~ is much larger than that of P~ and it is 
harder to steer a trajectory back into the basin P~ once 
it is in the basin of P^~ . 

Now we choose an initial condition such that the first 
few cells in the array are excited to the P^~ state (having 
a large basin of attraction compared to that of P~) and 
the rest are set to Pf state. In other words a wave front 
in the array is obtained by means of the two stable steady 
states. On actual numerical integration of Eqs.Q with 
N = 100 and with the diffusion coefficient chosen at a 
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FIG. 2: Space-time plot showing (a) propagation of wave 
fronts in one dimensional array (100 cells) of MLC circuits 
for D — 2.0 and (b) propagation failure for D — 0.4. 



higher value, D — 2.0, a motion of the wave front towards 
right (see Fig. |2Ia)) i s observed, that is a travelling wave 
front is found. After about 80 time units the wave front 
reaches the 100* n cell so that all the cells are now settled 
at the more stable state (P^) as demonstrated in Fig. 
121 a). When the value of D is decreased in steps and 
the analysis is repeated, the phenomenon of travelling 
wavefronts continues to be present. 

However, below a critical value of the diffusion coeffi- 
cient (D = D c ) a failure in the propagation has been 
observed, which in the present case turns out to be 
D = D c = 0.4. This means that the initiated wave- 
front is unable to move as time progresses and Fig. EJb) 
shows the propagation failure for D — 0.4. 



set of ten coupled first order odes: 





= j/i- h(xi) + 


D(x 2 
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= y 2 - h(x 2 ) + 
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X3 


= 1)2- h(x 3 ) + 
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+ X4), 


X4 


= 2/4 - h(xi) + 


D{x 3 


— 2X4 


+ x 5 ), 


X 5 


= 2/5 - h(x 5 ) + 


D{x 4 


- x 5 ), 




Vj 


= ~ u Vi ~ P x i, 


i = 


1,2.. 


,5, 



where h(x) is as given in Eq. and the parameters are 
fixed at /? = 1, a = 1, mo = —2.25, mi = 1.5, m 2 — 0.25 
and e = 0. 



1. Numerical analysis 

By carrying out a numerical analysis of the system 
(|4*|) with zero flux boundary conditions, we identify the 
following main results. 

a. We solve the above system of equations (0} for 
the chosen initial condition {xi(0), £2(0), x 3(0), 2:4(0), 
x 5 (0)} = {3.0, 3.0, -1.5, -1.5, -1.5} and ^(0) = -x,(0), 
i = 1, 2 ... 5, for different values of the diffusion coeffi- 
cient D in the region D £ (0, 1.0). Note that the values 
2:1(0) = 2; 2 (0) = 3.0 and 2:3(0) = 2:4(0) = 2:5(0) = -1.5 
correspond to the stable equilibrium points P^~ and P~ , 
respectively, of the uncoupled MLC oscillators. When D 
is decreased downwards from 1.0 we find that in the re- 
gion D > 0.4 propagation of wavefront occurs as in Fig. 
|2ta) (of the N = 100 cells system). At the critical value 
D = 0.4, the wavefront fails to propagate and asymptot- 
ically ends up in a stationary wavefront which is nearer 
to the initial condition (as in Fig. Hfb)). The wavefront 
failure phenomenon persists for all values of < D < 0.4. 

b. On the other hand, if we choose the initial con- 
dition as {2:1(0), 2:2(0), 2:3(0), 2:4(0), 2:5(0)} = {2.7990, 
2.1709, 1.7803, -1.4433, -1.4923} and 2/i(0) = -2^(0), 
i = 1,2, ..5, then the critical value D, where wavefront 
propagation failure occurs, turns out to be at a slightly 
higher value D = D c = 0.4703. 



2. Linear stability analysis 



B. Propagation failure mechanism: A case study of 
5 coupled oscillators 

To understand the failure mechanism, we have carried 
out a detailed case study of an array of five coupled MLC 
circuits, represented by the following set of equations. We 
have verified our conclusions for the case of 3, 4, 6 and 7 
oscillators also. The analysis can be in principle extended 
to arbitrary but finite number of oscillators. From Eq. 

, the present array can be represented by the following 



In order to understand the failure mechanism of wave- 
front propagation we consider all the possible steady 
states (equilibrium states) associated with Eqs.QJ for 
nonzero values of D. It is easy to check that Eq. Q 
possesses a maximum of 3 5 = 243 steady states. First 
we consider the three trivial steady states which can be 
easily obtained by considering x\ — x 2 — £3 = £4 = 2:5 
and yi = y 2 = y 3 = y 4 = 2/5, namely, X+ = 
(3.0, -3.0, 3.0, -3.0, 3.0, -3.0, 3.0, -3.0, 3.0, -3.0), 
X° = (0,0,0,0,0,0,0,0,0,0) and X~ 
(-1.5, 1.5, -1.5, 1.5, -1.5, 1.5, -1.5, 1.5, -1.5, 1.5). 
(Here the suffix p indicates that the steady states are 
the same as that of the D — case, namely P^~ , Pf and 
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From a linear stability analysis (details are given 
in Appendix B) we can easily conclude that the steady 
states X+ and X~ are stable while X® is unstable 
irrespective of the value of the diffusion coefficient D. 
Then any initial condition in the neighbourhood of Xjf 
will evolve into Xj^ asymptotically Also one has to note 
that, due to asymmetry in the function h(x), the basin 
of X+ is much larger than that of X~ . 

Since we are interested here to understand the propa- 
gation failure mechanism, we concentrate on those steady 
states which possess stationary wavefront structures near 
to the two types of chosen initial conditions in the nu- 
merical analysis discussed in Sec. IIII B II Out of the 
243 available steady states, we select a subset of them, 
X s = {X+, X°,X,r}, which is nearer to the numerically 
chosen initial state. These three steady states are defined 
in Ecis. ljBljl of Appendix B. They are identified such that 
the components x% s and X2 S are in the right extreme seg- 
ment of the characteristic curve, Fig. ^c), while xa s and 
X5 S are in the left extreme segment. On the other hand 
X3 S is allowed to be in any one of the three segments. 
The reason for such a choice is that we are looking for 
the formation of a wavefront which is closer to the cho- 
sen initial conditions in Sec. If 11 B II The corresponding 
forms of h(xi) as given in Eq. (J2J) are then introduced 
in Eqs.lQ} to find the allowed steady states. Then X+ 
corresponds to the case in which x^s as well as x\ s and 
X2s are in the right most segment while the rest are in 
the left most segment and X® and X~ correspond to the 
cases in which x^ s is, respectively, in the middle and the 
left extreme segments of the characteristic curve while 
X\si %2s and CC4 S , x$ s are chosen as in X+. These steady 
states may be explicitly obtained as a function of the dif- 
fusion coefficient D, whose forms are given in Appendix 
B. 

We now consider the linear stability of these steady 
states X+, X® and X~. The corresponding stability (Ja- 
cobian) matrices are obtained from the linearized version 
of Eq. , whose form is given in Eq. (|B4|> . The eigen- 
values are evaluated by numerical diagonalization of l|B4|) 
as a function D and the stability property analysed, ft 
has been found that out of these three steady states both 
X+ and X~ are stable for D < 0.4703 in which all of the 
eigenvalues of l|B4f) are having negative real parts, while 
X° is unstable for all values of D > 0. However, the 
real part of atleast one of the eigenvalues associated with 
each of X+ and X~ changes its sign to positive value 
at D = 0.4703 and thereby both the stationary fronts, 
X+ as well as X~, also lose their linear stability at this 
critical value. Thus for D > 0.4703, all the three steady 
states X+, X° and X~ are linearly unstable. 

First let us consider the second of the initial condi- 
tions chosen in Sec. If If B f I At and above the critical 
value D = 0.4703, all the three steady states are unsta- 
ble and so the system transits to the other nearby steady 
states. These are also found to be unstable by a similar 
analysis, ultimately ending up in the only nearby avail- 
able stable state which is X+ = ((3.0, -3.0), (3.0, -3.0), 



(3.0,-3.0), (3.0,-3.0), (3.0,-3.0)), thereby explaining 
the propagation of the wavefront starting from the initial 
condition X(0) = {(2.7990, -2.7990), (2.1709, -2. f 709), 
(1.7803,-1.7803), (-1.4433,1.4433), (-1.4923,1.4923)} 
(which is the second of the numerical results given in 
Sec. HUB 111 . However for D < 0.4703, since X+ is sta- 
ble and is also quite closer to the initial state, the system 
settles down in this state itself and no propagation oc- 
curs, thereby explaining the phenomenon of propagation 
failure as a subcritical bifurcation when the diffusion co- 
efficient D is increased from upwards. Note that even- 
though X~ is also stable, it is far away from the initial 
state and so no transition to this state will occur. A ta- 
ble of these steady states X+ and X~ as a function of D 
(Table ^| is given confirming that X+ is always nearer 
to the initial state compared to X~ . 

Now let us consider the first of the initial conditions 
of our numerical analysis (see Sec. IIII B 111 , where wave- 
front propagation occurs at D = 0.4 itself and not at 
D = 0.4703, eventhough the steady states X+ and X~ 
are still linearly stable in the region 0.4 < D < 0.4703. 
In order to understand this aspect, we start analysing 
the nature of the basin of attraction associated with 
the steady state X~ for different D values numerically 
(Note that for the present set of initial conditions X~ 
is nearer to it than X+). We find that the basin of 
X~ shrinks as D is increased from upwards and van- 
ishes completely at D — 0.4703. We also find that the 
chosen initial condition, X(0) = ((3.0, -3.0), (3.0, -3.0), 
(-1.5, 1.5), (-1.5, 1.5), (-1.5, 1.5) falls fully within the 
basin of X~ as long as D < 0.4. However due to the 
shrinking nature of the basin as D increases (as given 
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FIG. 3: Figure showing partial basin for the steady state 
X~ for various values of the diffusive coupling coefficient. In 
the figure each line corresponds to an initial condition. The 
figure is drawn by fixing the initial state of the first oscillator 
at xi (0) = 3.0, yi(Q) = —3.0 and for clarity decreasing the 
rest downwards in steps of 0.2. 

in Figs. ©), a part of the initial state falls outside the 



6 



basin for D > 0.4 (Figs. Etb)-(d)). As a result the cho- 
sen initial condition does not end up in the stationary 
front X~ for 0.4 < D < 0.4703 indicating a global in- 
stability (eventhough the state is linearly stable) . There- 
fore propagation starts to occur as soon as D > 0.4, 
while it fails for D < 0.4 and it finally ends in the 
stable state A+ described in the previous case. How- 
ever, if an initial condition which happens to fall com- 
pletely within the basin of attraction of X^ 1 , the prop- 
agation starts to occur only for values of D > 0.4703. 
In fact this is what happens for the second chosen initial 
condition (X(0) = {(2.7990, -2.7990), (2.1709, -2.1709), 
(1.7803, -1.7803), (-1.4433, 1.4433), (-1.4923, 1.4923)}, 
as explained above (Note that the basin structure shown 
in Fig. does not cover this case). 
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FIG. 4: Three dimensional plot showing propagation of two 
dimensional wave fronts (100 x 100) at various time units (a) 
T = 0, (b) T = 50, (c) T = 100 and (d) T = 125 for Di = 2.0 
and Di = and (e) the same plot for propagation failure case 
Di = 0.4 and D 2 = at T = 1000 t.u. 

In conclusion, the propagation failure of wave fronts in 
arrays of coupled nonlinear oscillators is essentially due 
to the existence of stationary fronts (stable steady states) 
near to the initial states for low diffusion coefficients, be- 
low a certain critical value, which lose stability via sub- 
critical bifurcation. This coupled with the existence of 
necessary basin of attraction (for a global stability) of 
the steady state gives rise to the specific critical value 
of D = D c above which propagation begins for a given 
initial condition. 



C. Propagation phenomenon in two-dimensional 
array 

For studying the propagation phenomenon in two di- 
mensional array of coupled MLC circuits, we have con- 
sidered the case in which F — 0, D\ ^ and D2 = 0. 
The propagation of wave front and its failure have been 
observed as in the case of one dimensional array. A two 
dimensional array of 100 x 100 coupled MLC circuits has 
been considered for the numerical simulation with initial 
conditions chosen such that the first few cells in the array 
are excited to the state and the other cells are set at 
the P~ state. Figs. EIa)-0Jd) show the propagation of 
the wave front in the two dimensional array for D\ = 2.0 
and D2 = at various time units. Again, below a cer- 
tain critical diffusion coefficient, D\ = 0.4, propagation 
failure occurs. One such case is illustrated in Fig. 0{e) 
for Di = 0.4 and D% = 0. The phenomenon can again 
be explained through a linear stability analysis as in the 
case of the one dimensional array. 

D. Effect of weak coupling 

In the above, the investigation has been made by con- 
sidering the system as an ideal one (as far as the circuit 
parameters are concerned). But from a practical point of 
view, there are defects in the coupling parameters which 
may result in a weak coupling at any of the cell in the 
array. In the following we study the effect of such weak 
coupling on the propagation of wave front. 

Let us consider a weak coupling at the fc^h cell. By 
this we mean that the k l cell in the array is coupled 
to its nearest neighbour (k + l)" 1 and (k — l)^ n cells 
by resistors with slightly higher values than that of the 
others. We have studied the effect of this defect on the 
propagation of wave fronts. Numerical simulations have 
been carried out by considering an one dimensional ar- 
ray with 100 cells, where the initial conditions are chosen 
as in the case of propagation phenomenon in regular one 
dimensional arrays (see Sec. IIII Af) . From the numerical 
simulation results, we find that there is an abrupt stop in 
the propagation when the wave front reaches the weakly 
coupled cell. This happens when the coupling coefficient 
on either side of the cell in the array has a value even 
above the critical value (D — 0.4) for propagation failure 
discussed in Sec. IIII Al Fig. [5] shows a blocking in the 
propagation of the wave front when the coupling coeffi- 
cients on either side of the fc^ n cell (k = 25), which we 
call as Dk, is set to 0.47173 with the rest of the coupling 
coefficients set to 1 (D = 1). We observe that the ac- 
tual blocking occurs when the wave front reaches the k^ 1 
cell. We can say this is a kind of failure in the propaga- 
tion (because the wavefront will never reach the last cell 
(that is the lOO^h cell in the array) by means of blocking. 

As far as the initial conditions are concerned, since the 
system is in the propagation region (D = 1.0) any choice 
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cell will block the propagating front and hence it will 
end up in the stationary front X~ asymptotically for 
D k < 0.4717. For D k > 0.4717, the stationary state 
X~ is no longer stable and hence the propagation occurs 
without any blocking. Thus the blocking of the wave 
front can also be attributed to the existence of a sta- 
ble stationary front when weak coupling is present in the 
system. 



i 



FIG. 5: Figure showing the effect of weak coupling at the 
25 th cell in the one dimensional array (Eq. (3)). (a) three 
dimensional space-time plot and (b) density plot. 

of initial condition of wavefront structure will lead to uni- 
form propagation upto the weakly coupled cell. Now our 
motivation is to find the property of the above uniform 
wavefront when it reaches the weakly coupled cell. 

One can analyse the above phenomenon by means of a 
linear stability analysis of the steady states as in the case 
of propagation phenomenon in Sec. IIII Bl We consider 
again as an example the system of five coupled MLC 
circuits as in Eqs.QJ in which the middle cell (x^) is 
now taken as an example to be coupled weakly to its 
neighbours, respectively x 2 and x±. This set up can be 
represented by the following set of ten coupled odes: 



X\ 


= V\ - h(xi) + 


D(x 2 - 


Xl), 






±2 


= Vi - h{x 2 ) + 


D(xi - 


xa) - 


\-D k (x 3 


- x 2 ), 


is 


= y 2 - h(x 3 ) + 


Dk{x2 


-2x 3 


+ Xi), 




i'4 


= 2/4 - h(x 4 ) + 


Dk(x 3 


- x 4 ) 


+ D(x 5 


- x 4 ), 


x 5 


= y 5 - h(x 5 ) + 


D[xa - 


X5), 






>jj 


= -oVj ~ faj, 


J = l 


,2... 


,5, 





Here D k corresponds to the coupling coefficient of the 
weakly coupled cell and other parameters are the same 
as in QJ. As in Sec. IIII Bl we calculated the steady 
states associated with (JjJ and analysed the linear stabil- 
ity properties. 

Let us now consider the three steady states of the sys- 
tem (0), namely, X s — {X+,X®,X~} which possess 
wavefront structures near the chosen initial conditions. 
These steady states can be found as discussed in Sec. 
IIII Bl By fixing D = 1.0, we analyse the stability of X s 
as a function of D k . We find that, the steady states, Xf- 
are linearly stable for values of Dk, < Dk < 0.4717 
and become unstable for D k > 0.4717 via subcritical bi- 
furcation. However X® state continues to be unstable 
irrespective of the value of Dk ■ 

At this point, one may think of an array with large 
number of MLC circuits. Since all the cells, except the 
k^ 1 cell, have diffusion coefficient D = 1.0 which is in the 
propagation region, the appropriately chosen initial con- 
dition will evolve as a propagating front upto the weakly 
coupled cell. When this wave front reaches the weakly 
coupled cell, the existence of the stationary front A s 7 
for the array due to the presence of the weakly coupled 



E. Turing patterns 

Another interesting dynamical phenomenon in the cou- 
pled arrays is the formation of Turing patterns. These 
patterns are observed in many reaction diffusion systems 
when a homogeneous steady state which is stable due 
to small spatial perturbations in the absence of diffu- 
sion becomes unstable in the presence of diffusion[Turing, 
1952]. To be specific, the Turing patterns can be ob- 
served in a two variable reaction-diffusion system when 
one of the variables diffuses faster than the other and 
undergo Turing bifurcation, that is, diffusion driven in- 
stability [Turing, 1952; Murray, 1989]. 

Treating the coupled array of MLC circuits as a dis- 
crete version of a reaction-diffusion system, one can as 
well observe the Turing patterns in this model also. For 
this purpose, one has to study the linear stability of sys- 
tem (3) near the steady state. In continuous systems, the 
linear stability analysis is necessary to arrive at the con- 
ditions for diffusion driven instability. A detailed deriva- 
tion of the general conditions for the diffusion driven 
instability can be found in Murray [1989]. For discrete 
cases one can follow the same derivation as in the case of 
continuous systems by considering solutions of the form 
expi(kj — Xt) [Murray 1989; Munuzuri et ai, 1995]. Here 
k and A are considered to be independent of the position 
j (j = 1,2, ••• ,N). For Eq. (3), the criteria for the 
diffusion driven instability can be derived by finding the 
conditions for which the steady states in the absence of 
diffusive coupling are linearly stable and become unstable 
when the coupling is present. One can easily derive from 
a straightforward calculation that the following condi- 
tions should be satisfied for the general reaction-diffusion 
system of the form given by Eq. 

fx + 9y < 0, 
Ix9y ~ f V 9x > 0, 

f x D 2 -g y D 1 > 0, 
(f x D 2 - g v D{f - ±D x D 2 (f x D 2 - g y D x ) > 0. (6) 

The critical wave number for the discrete system (3) can 
be obtained as 

Combining Eqs. ©-0, one obtains [Munuzuri et al., 
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1995] the condition for the Turing instability such that 
fxD2 - g y Di 



8£>i£>2 



< 1. 



(8) 



We apply these conditions to the coupled oscillator 
system of the present study. For this purpose, we fix 
the parameters for the two-dimensional model (3) as 
{/3,cr,e,TOo,mi,m 2 } = {0.8, 0.92,0.1,-0.5,0.5,0.5} and 
verify that this choice satisfies the conditions iJBJ to JSJ- 
The numerical simulations are carried out using an ar- 
ray of size 100 x 100 and random initial conditions near 
the steady states are chosen for the x and y variables. 






FIG. 6: The spontaneous formation of Turing patterns in 
an array of 100 x 100 oscillators for the parameters f3 = 0.8, 
a = 1, m = -0.5, mi = 0.5, m 2 = 0.5, e = 0.1, F = 0.0, 
uj = 0.75, Di = 1 and D2 = 10 in Eq. (3) at various time 
units (a) T = 50, (b) T = 100, (c) T = 500 and (d) T = 2000 
(e) the same figure (d) in the Fourier space. 

Figs. EJa)-UJd) show how the diffusion driven instabil- 
ity leads to stable hexagonal pattern (Fig. Efd)) after 
passing through intermediate stages (Figs. Et a )"Et c ))- 
Fig. Ele) shows the two dimensional Fourier spectrum 
of the hexagonal pattern (Fig. Eld)). Further, the spon- 
taneously formed patterns are fairly uniform hexagonal 
patterns having a penta-hepta defect pair. These defects 
are inherent and very stable. We note that such Turing 
patterns have already been observed in an array of cou- 
pled Chua's circuits [Munuzuri et al, 1995]. However, our 



aim here is to investigate the effect of external forcing on 
these patterns, which is taken up in the next section. 



IV. SPATIO-TEMPORAL PATTERNS IN THE 
PRESENCE OF PERIODIC EXTERNAL FORCE 

The effect of external fields on a variety of dynamical 
systems has been studied for a long time as driven sys- 
tems are very common from a practical point of view. For 
example, in a large number of dynamical systems includ- 
ing the Duffing oscillator, van der Pol oscillator and the 
presently studied MLC circuits, temporal forcing leads 
to a variety of complex dynamical phenomena such as 
bifurcations, quasiperiodicity, intermittency, chaos and 
so on[Guckcnhcimcr & Holmes, 1983; Hao Bai-Lin, 1984; 
Lakshmanan & Murali,1996]. Also the studies on the ef- 
fect of external fields in spatially extended systems have 
been receiving considerably increasing interest in recent 
times [Pismen, 1987; Walgraef, 1988; 1996]. Particularly, 
with the recent advances in identifying localized and os- 
cillating structures and other spatio-temporal patterns 
in driven nonlinear dissipative systems such as granular 
media, driven Ginzburg-Landau equations and so on, it 
is of considerable interest to study the effects of forcing 
on array of coupled systems such as (3). Motivated by 
the above, we investigate the effect of external forcing 
in the propagation of wave front and formation of Tur- 
ing patterns in the coupled MLC circuits in one and two 
dimensions. 



Effect of external forcing on the propagation of 
wave fronts 



In this subsection we study the effect of external forc- 
ing on the propagation of wave fronts. For this purpose 
we consider an one dimensional array of coupled MLC cir- 
cuits with initial and boundary conditions as discussed 
in Sec. IIII Al Now we perform the numerical integration 
by the inclusion of external periodic force of frequency 
u> = 0.75 in each cell of the array (see Eq. (1)). By vary- 
ing the strength, F, of the external force we study the 
behaviour of the propagating wave front in comparison 
with the force free case (F = 0) as discussed in Sec. 1111 Al 
We find that in the propagation region (D > D c ), the ef- 
fect of forcing is just to introduce temporal oscillations 
and the propagation continues without 

any disturbance (see Fig. Efa)) as in the case F = 
(Fig. |2[a)). Of course this can be expected as the exter- 
nal force is periodic in time. However, interesting things 
happen in the propagation failure region discussed in Sec. 
IIII Al In this region, beyond a certain critical strength of 
the external forcing, the wavefront tries to move a little 
distance and then stops, leading to a partial propagation. 
Fig. |7f b) shows such a partial propagation observed for 
F = 0.6 and D = 0.22 (This may be compared with Fig. 
I21b)). The phenomenon can be explained by considering 
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(a) <b| 




FIG. 7: Propagation phenomena in the presence of forcing: 

(a) Propagation of wavefront for D = 2.0 and F = 0.6 and 

(b) The partial propagation observed for the D — 0.22 and 
F = 0.6. 

the propagation failure mechanism discussed in Sec. IIII Bl 
in which one may look for a spatially stationary and tem- 
porally oscillating wavefront. The initial wavefront tries 
to settle in the nearby stationary state. However, the 
system will take a little time and space to settle due the 
effect of forcing combined with the transient behaviour 
of the system. Thus the inclusion of external forcing in 
the propagation failure region can induce the wavefront 
to achieve a partial propagation. 



B. Transition from hexagons to rhombs 

It is well known that the defects are inherent in very 
many natural pattern forming systems. In most of the 
pattern forming systems, the observed patterns are not 
ideal. For example, the patterns are not of perfect rolls 
or hexagons or rhombs. A commonly observed defect in 
such systems is the so called penta-hepta defect (PHD) 
pair which is the bound state of two dislocations [Tsim- 
ring, 1996]. Experiments on spatially extended systems 
often show the occurrence of PHD in spontaneously de- 
veloped hexagonal patterns [Pantaloni & Cerisier, 1983]. 
In the present case also, the existence of PHD pair can 
be clearly seen from Fig. HJd). In such a situation, it is 
important to study the effect of external periodic force 
in the coupled arrays of MLC circuits. 

Now we include a periodic force with frequency w and 
amplitude F in each cell of the array and we numeri- 
cally integrate Eqs.Q using fourth order Runge-Kutta 
method with zero flux at boundaries. By fixing the fre- 
quency of the the external periodic force asw = 0.75 and 
varying the amplitude (F) we analyse the pattern which 
emerges spontaneously. Interestingly for F = 0.25, the 
defects (PHD pair) which are present in the absence of 
external force (Fig. EJd)), gets removed resulting in the 

transition to a regular rhombic array. Fig. [8(a) 
shows the gray scale plot of the pattern observed for 
F = 0.25 and the corresponding plot in Fourier domain 
(Fig. H[b)). Thus, from the above we infer that the in- 
clusion of external periodic force can cause a transition 
from hexagonal pattern to rhombic structures. We note 
that Perez-Munuzuri et al. have observed, in arrays of 
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FIG. 8: Figure showing the formation of rhombic structures 
in the presence of forcing, (a) Gray scale plot of the rhombic 
pattern for F — 0.25 and (b) The same figure in Fourier space. 

coupled Chua's circuits, similar transition from hexagons 
to rhombs by means of sidewall forcing [Perez-Munuzuri 
et al., 1995]. However in our case we apply forcing simul- 
taneously to all the cells, in order to mimic situations 
such as application in vertically vibrated granular me- 
dia [Umbanhowar et al., 1996]. 

C. Transition from hexagons to rolls 

In addition to the transition from hexagons to 
rhombs by the influence of external periodic force, there 
are also other possible effects due to it. To real- 
ize them, we consider a different set of parametric 
choice {P, a, m , mi, m 2 } = {0.734722, 0.734722, -0.874, 
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FIG. 9: Figure showing the transition from hexagons to rolls 
for {/3, a, e, m , mi, m 2 } = {0.734722, 0.734722, 0.15, -0.874, 
-0.4715, -0.4715} with Di = 1 and D 2 = 5. (a) F = 0, (b) 
F = 0.15, (c) F = 0.35 and (d) F = 0.45. 

-0.4715, -0.4715} with e = 0.15, D x = 1.0 and D 2 = 
5.0. For this choice the system shows hexagonal patterns 
with defects including domains of small roll structures 
(Fig. Ha)). 

Now when the external periodic force is included a 
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FIG. 10: Snapshots showing the breathing oscillations for 
F = 0.05 and {/?, a, e, mo, mi, m 2 } = {0.734722, 0.734722, 
0.15, -0.874, -0.4715, -0.4715} with £>i = 2 and D 2 = 5 for 
various time units starting from T = 4001. 



0.734722, 0.10, -0.874, -0.4715, -0.4715} with D x = 2.0 
and £>2 = 5.0 such that a regular hexagonal pattern is 
observed in the absence of external periodic force. From 
numerical simulations, we observed that a space-time pe- 
riodic oscillatory pattern (breathing motion) sets in for 
a range of low values of F. Fig. ^| shows the typical 
snapshots of the oscillating pattern at various instants 
for the specific choice F — 0.05. We have integrated over 
10000 time units and the figure corresponds to the re- 
gion T = 4000 - 4018. Typically we find that the breath- 
ing pattern repeats itself approximately after a period 
T = 15.0 in the range of our integration. One may con- 
clude that the emergence of such breathing oscillations 
is due to the competition between the Turing and Hopf 
modes in the presence of external periodic force. 



V. SPATIO-TEMPORAL CHAOS 

Next we move on to a study of the spatio-temporal 
chaotic dynamics of the array of coupled MLC circuits 
when individual cells are driven by external periodic 
force. The motivation is that over a large domain of 
(F, oj) values the individual MLC circuits typically ex- 
hibits various bifurcations and transition to chaotic mo- 
tion (see Appendix A). So one would like to know how 
the coupled array behaves collectively in such a situa- 
tion, for fixed values of the parameters. For this pur- 
pose, we set the parameters at {/3, a, e, too, mi, TO2, w} = 
{1.0, 1.015, 0, -1.02, -0.55, -0.55, 0.75}. In our numeri- 
cal simulations, we have mainly considered the one di- 
mensional array specified by Eq. (1) and assumed peri- 
odic boundary conditions. The choice of periodic bound- 
ary conditions here makes the calculation of Lyapunov 
exponents easier so as to understand the spatio-temporal 
chaotic dynamics of (1) in a better way. 



transition in the pattern from hexagonal structure to rolls 
starts appearing. By fixing the frequency of the external 
periodic force again at oj = 0.75, we observed the actual 
transition from hexagons to rolls as we increase the forc- 
ing amplitude (F). Figs. Etb)-|H{d) show the gray scale 
plots for F = 0.15, 0.35 and 0.45, respectively. Obvi- 
ously the transition is due to the existence of small roll 
structures in the pattern for F = which nucleates the 
formation of rolls in the presence of forcing. 

D. Breathing oscillations 

In the above, we have shown that the inclusion of the 
external periodic force can make a transition from one 
stationary pattern to another stationary pattern like the 
transition from hexagons to rolls. Besides these, are there 
any time varying patterns? As mentioned above, pat- 
terns such as localized and breathing oscillations have 
considerable physical interest. In this regard, we consid- 
ered the parameters {/3, a, e, mo, mi, = {0.734722, 



A. Spatio-temporal regular and chaotic motion 

Numerical simulations were performed by consider- 
ing random initial conditions using fourth order Runge- 
Kutta method for seven choices of F values. The cou- 
pling coefficient in Eq. (1) was chosen as D = 1.0. Out of 
these, the first three lead to period-T, period-2T, period- 
4T oscillations, respectively and the remaining choices 
correspond to chaotic dynamics of the single MLC cir- 
cuit. Figs. [TlTal- triT g') show the space-time plots for 
F = 0.05, F = 0.08, F = 0.09, F = 0.12, F = 0.13, 
F = 0.14 and F = 0.15, respectively. From the Figs. 
ETTal-OTTo). it can be observed that for F = 0.05, 0.08 
and 0.09 the MLC array also exhibits regular periodic be- 
haviour with periods T, 2T and 4T respectively in time 
alone as in the case of the single MLC circuit. However, 
for F = 0.12 and 0.13(Figs. UTTd) andCJc)) one ob- 
tains space-time periodic oscillations eventhough each of 
the individual uncoupled MLC circuits for the same pa- 
rameters exhibits chaotic dynamics. We may say that a 
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FIG. 11: Spatio-temporal periodic oscillations in a grid of 
100 oscillators for the external periodic forcing strength (a) 
F = 0.05, (b) F = 0.08, (c) F = 0.09, (d) F = 0.12, (e) 
F = 0.13, (f) F = 0.14 and (g) F = 0.15 



B. Size dependence of the ST Chaos 

Since the above study of spatio-temporal chaos in- 
volves a large number of coupled chaotic oscillators it 
is of great interest to analyse the size dependence of the 
dynamics of these systems. To start with, we consider 
the case of 10 coupled oscillators with periodic boundary 
conditions and numerically solve the system using fourth 
order Runge-Kutta method with the other parameters 
chosen as in Sec. IV Al The value of F is chosen in the 
range (0.12,0.15). We find that this set up shows differ- 
ent behaviour as compared to the 100 cell case. Actually 
the system gets synchronized to a chaotic orbit rather 
than showing periodic behaviour as in the case of 100 
cells described above. 

First we analyse the dynamics for F = 0.12 by slowly 
increasing the system size from N = 10. We noted 
that for the system size, N < 42 it gets synchronized 
to a chaotic orbit as mentioned above. We have also 
calculated the corresponding Lyapunov spectrum and 
found that there is one positive maximal Lyapunov ex- 
ponent as long as N < 42. For example, for N = 42, 
Amax = 0.1162, with the rest of the exponents being 
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FIG. 12: (a) The chaotic attractor at the 5 L11 cell for the 
synchronized state (N = 42) and (b) the periodic orbit in the 

5 th cell for the controlled state (N = 43) 



kind of controlling of chaos occurs due to the coupling, 
though the coupling strength is large here. From the 
above analysis it can be seen that the macroscopic sys- 
tem shows regular behaviour in spite of the fact that 
the microscopic subsystems oscillate chaotically. Finally 
for F — 0.15 the coupled system shows spatio-temporal 
chaotic dynamics (Fig. Illf g)) and this was confirmed by 
calculating the Lyapunov exponents using the algorithm 
given by Wolf et al. [Wolf et al., 1985]. For example, we 
calculated the Lyapunov exponents for N = 50 coupled 
oscillators (For N = 100, it takes enormous computer 
time relatively which we could not afford at present) 
and we find the highest three exponents have the val- 
ues Amax = A = 0.1001, Ai = 0.0776, A 2 = 0.0092 
and the rest are negative (see next subsection for further 
analysis). We also note that for F = 0.14, the system is 
in a transitional state as seen in Fig. Illf f). 



negative. Fig. I12f a) shows the dynamics of the 5 cell 
in the array and !13f a) depicts the space-time plot of the 
synchronized spatio-temporal chaos. The system shows 
entirely different behaviour when we increase the system 
size to N = 43. As noted in the previous subsection, 
IV Al there occurs kind of suppression of spatio-temporal 
chaos. Fig. I12f b) shows the resultant periodic orbit in 
the 5^ n cell of the array and Fig. I13f b) shows the space- 
time plot of the spatio-temporal periodic behaviour of the 
array. The maximal Lyapunov exponent is found to be 
negative in this case (Amax = —0.001474). We observed 
the same phenomenon for F = 0.13 also. 

Next we consider the case of F = 0.15 in which 
the coupled system in the previous subsection showed 
spatio-temporal chaos. From numerical simulations, we 
again observed a synchronized motion for N < 31 and 
the corresponding Lyapunov spectrum shows one posi- 
tive exponent only with all other exponents being nega- 
tive. But for N > 31, the coupled system shows spatio- 
temporal chaos. The Lyapunov spectrum in this case (for 
N > 31) possesses multiple positive exponents. For ex- 



12 



1.12 
1 118 
,1.116 
*t.1M 
1.112 
1 11 




R i L L 

i — M/^-^m\ 



FIG. 13: Space-time plot showing the size dependence of the 
dynamics in the coupled chaotic oscillators, (a) N = 42 and 
(b) N = 43. 



ample, for N = 43, (A ma x = A = 0.0997, A x = 0.0633, 
A 2 = 0.0038). 

From the above analysis we sec definite evidence that 
the spatio-temporal chaotic dynamics in coupled MLC 
circuits depends very much on the size of the system. 
A similar analysis can be performed to investigate the 
nature of the spatio-temporal dynamics on the strength 
of the coupling coefficient D. Preliminary analysis shows 
similar phenomenon at critical values of D. Further work 
is in progress on these lines. 



VI. CONCLUSIONS 

In this paper we have discussed various spatiotemporal 
behaviours of the coupled array of Murali-Lakshmanan- 
Chua circuits in one and two dimensions. In the first 
part we have reported the propagation phenomenon of 
the travelling wave fronts and failure mechanism in the 
absence of external forcing. We have shown that the 
propagation phenomenon is due to the loss of stability 
of the steady states via subcritical bifurcation coupled 
with the presence of necessary basin of attraction of the 
steady states. The effect of weak coupling on the travel- 
ling wave phenomenon has also been studied. In the two 
dimensional array we studied the onset of Turing insta- 
bility leading to the spontaneous formation of hexagonal 
patterns. 

Introduction of external periodic force has been shown 
to lead to a partial propagation in the failure region of one 
dimensional array, while in the two dimensional array a 
transition from hexagons to rhombs take place. We have 
also showed that a transition from hexagons to rolls can 
occur by the influence of external periodic force provided 
domains of small roll structures already present in the ab- 
sence of external force. Space-time breathing oscillations 
as a result of the competition between Hopf and Turing 
modes in the presence of external periodic force have also 
been observed. We have also studied the spatio-temporal 
chaotic dynamics and showed that the MLC array ex- 
hibits both periodic and chaotic dynamics as in the case 
of single oscillator. In the chaotic regime the dynamics is 
also critically dependent on the system size. For a large 
array, the collective behaviour shows periodic oscillations 
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FIG. 14: A single MLC circuit 



both in space and time eventhough the individual sub- 
systems oscillate chaotically as the coupling brings the 
macroscopic system into a regular motion. On the other 
hand, below a critical size of the system, it gets synchro- 
nized to a chaotic orbit. The study of array of diffusively 
coupled driven nonlinear oscillators is of intrinsic interest 
as they represent many natural phenomena such as Fara- 
day instability, patterns in granular media and so on. Of 
particular interest will be to look for localized structures 
in these systems. Whether such excitations exists in the 
present arrays of oscillators is a question to be answered 
in the near future. Work is in progress along these lines. 
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APPENDIX A: DYNAMICS OF THE MLC 
CIRCUIT 



The Murali-Lakshmanan-Chua circuit, Fig. I14f a). is 
the simplest second order dissipative nonautonomous cir- 
cuit, consisting of Chua's diode as the only nonlinear 
element [Murali et al, 1994a]. This circuit contains a ca- 
pacitor (C), an inductor (L), a linear resistor (R), an ex- 
ternal periodic forcing (/ sin fit) and a Chua's diode. By 
applying the Kirchoff's laws to this circuit, the governing 
equations for the voltage v across the capacitor C and the 
current i l through the inductor L are represented by the 
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following set of two first order nonautonomous differen- 
tial equations: 



c Tt = iL ~ 9{v) > 

L^- = -Ri L - R s i L - v + /sin fit, (Al) 

where g(v) is a piecewise linear function corresponding 
to the characteristic of the Chua's diode (N) and is given 
by 

f e' + G 2 v R + (G - Gi), v R > B p 
9( v r) = { e ' + Govr, -B p <v R <B p 

(e' + dvR-iGo-Gi), vr<-B p 

The piecewise nature of the characteristic curve of 
Chua's diode is obvious from Eq. I|A2(1 [Murali et al, 
1994a]. The slopes of left, middle and right segments of 
the characteristic curve are Gi , Go and G2 , respectively. 
—Bp and B p are the break points and e' corresponds to 
the dc offset in the Chua's diode. Rescaling Eq. (|A1|) as 
v = xB p , i L = GyBp, G = 1/R,u> = fiG/G t = rC/G 
and e = e'/G and then redefining t as t the following set 
of normalized equations are obtained: 

x = y- h{x), 

y = — fix — ay + F sin uit, (A3) 

with 

{e + m 2 x + (to — Tfi\ ) , x > x 2 
e + mox, x\ < x < x 2 , (A4) 

e + m\X — (mo — mi), x < X\ 

where fi = (C/LG 2 ), a = (C/LG 2 )(l + GR S ) and F = 
(ffi/Bp). Obviously h(x) takes the form as in Eq. 
with too = Go/G, mi = Gi/G and m 2 = G 2 /G. The 
dynamics of Eq. (|A3(I depends on the parameters fi, a, 
too, toi, m 2 , e, lo and F. 

The rescaled parameters corresponding to the experi- 
mental observations reported in Murali et al. [1994a] cor- 
respond to fi = 1.0, a = 1.015, too = —1.02, m\ = —0.55, 
to 2 — 0.55, e = and u> = 0.75. By varying F one can 
observe the familiar period doubling bifurcations lead- 
ing to chaos and several periodic windows in the MLC 
circuit. Fig. I15f a) shows the one parameter bifurcation 
diagram in the F-x plane [F € (0,0.7)]. A summary of 
the bifurcations that occur in this case for different F val- 
ues is given in Table HI Further it is of great interest to 
consider the parametric choice {too, mi, m 2 , e, fi, a, lu} = 
{-2.25, 1.5, 0.25, 0.0, 1.0, 1.0, 0.75} which corresponds to 
the function h(x) having the form as shown in Fig. 
dc). This choice of parameters provides the possibil- 
ity of bistability nature in the asymmetric case in the 
absence of periodic forcing. In this case one can easily 
observe from numerical simulations that the MLC circuit 
admits only limit cycles for F E (0, 0.7). The bifurcation 
diagram in the F-x plane is depicted in Fig. I15f b). 



FIG. 15: Bifurcation diagram in the F-x plane of the MLC 
circuit for the parametric choice (a) {fi, a, mo, mi, mi, e, u>} 
= {1.0,1.015,-1.02,-0.55,-0.55,0,0.75} and (b) {fi,a,m , 
mi, m 2 , e, u} = {1.0, 1.0, -2.25, 1.5, 0.25, 0, 0.75}. 



TABLE I: Summary of bifurcation phenomena of Eq.(A3) 
amplitude (F) description of attractor 
< F < 0.071 period-1 limit cycle 

0.071 < F < 0.089 period-2 limit cycle 
0.089 < F < 0.093 period-4 limit cycle 
0.093 < F < 0.19 chaos 
0.19 < F < 0.3425 period-3 window 
0.3425 < F < 0.499 chaos 
0.499 < F < 0.625 period-3 window 
0.625 < F period-1 boundary 

In the main text, we also consider other choices of the 
parameters depending on the phenomenon we are inter- 
ested in. The corresponding values are given at the ap- 
propriate places in the text. 



APPENDIX B: STEADY STATES AND 
JACOBIAN MATRIX OF THE SYSTEM OF 5 
COUPLED OSCILLATORS (EQS. ©) 

The steady states of Eq. (@J can be obtained ex- 
plicitly as a function of the diffusion coupling coeffi- 
cient D by introducing the appropriate form of h(xi) 
from Eq. J3J. Out of the 3 5 = 243 possible steady 
states we consider here a subset X s = {X+ , X ( J , X~} 
and X p = {X+, X°, X- } discussed in Sec. HUB 21 

To obtain the set X s we proceed as follows. The com- 
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ponents x\ s and Xi s are in the right extreme segment of 
the characteristic curve, Fig. ^c), while x^ s and x§ s are 
in the left extreme segment. Therefore choosing h(xi) = 
m2X + {mo~mi) for i = 1,2 and h(xi) — mix — {mQ — mi) 
for i = 4, 5 (see Eq. we obtain 



'4s 



"5s 



-^(64L> 4 + 1280D 3 + 350UD 2 + 2750L> + 625), 
-^(64L> 4 + 1760Z) 3 + 3800L> 2 + 2750L> + 625), 



1,2,. 



.5, 



-^"s" — i x ts> Vlsi X 2s> Utsi X ts>yts' X ts' Vis' X ts'V5s} 



(Bla) 



for h(xs) = mix + (mo — m-i), 

X-a = i x ls j 2/ls ) x 2s > 2/2s ) > V3s > ^4^ ) 2/4s > x 5s > Vhs } 

for h{xz) = vtiqx and 



(Bib) 



-^s — l^ls' 2/ls' X 2s' 2/2s' ^s' 2/3s' ^48' 2/4 s' %si 

for /i(a;3) = mix— (rrio — m i)- Explicitly solving Eqs. Q 
we can write 



— (16D 3 + 620Z? 2 + 55CLD + 125), 
K 



'2s 



— (16L> 3 + 500D 2 + 550D + 125), 
K 

4=(16D 3 + 260L> 2 + 400L> + 125), 
K 

— (32D 3 - 200L» 2 - 400L> - 125), 
K 

^—(32L> 3 - 560L> 2 - 550L> - 125), 
K 

— X; , i = 1, 2, ... 5, 



with 



x 2s 



1 3s 



1 4 s 



' 5 s 



15 

— (64L> 3 + AAD 2 



50D - 21 
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M 
30D 



^■(48L> 3 + AAD 2 - SOD - 25), 



(8£> 2 + 12D + 5), 



- — (32D 3 - 32D 2 - 40L> - 25), 

_Z^(64D 3 + 32L> 2 - 50D - 25), 

-xl i = l,2,...5. (B3) 



K = 112D 3 + 620D 2 + 550£> + 125. 



L = 256L> 4 + 1880L> 3 + 3800L> 2 + 2750Z) + 625. 



'-2s 



! 3s 



--(32D 4 - 920L» 3 - 3200L> 2 - 2750Z) - 625), 
L 

-y(32D i - 560D 3 - 2000L> 2 - 2000Z) - 625), 
L 

- — (64D 4 + 320L> 3 + 1700L> 2 + 2000Z? + 625), 

(B2) 



and 



M = 64L> 4 + 320L> 3 + 140L> 2 - 250L> - 125. 



The stability determining eigenvalues are the roots of 
the characteristic equations of the Jacobian matrix for 
the linearised version of Eq. Q: 



f-h 



J{X S ) 





D 











1 











o\ 


D 


-h X2 - 2D 


D 











1 














D 


-h X3 - 2D 


D 











1 














D 


-h X4 - 2D 


D 











1 














D 


-h X5 - D 














1 


-1 














-1 

















-1 














-1 

















-1 














-1 

















-1 














-1 

















-1 














-l) 



(B4) 



where h x is the derivative of h(x) with respect to x eval- uated at the steady states. 
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Particularly we can show that for the trivial steady 
states Xf = (3.0, -3.0, 3.0, -3.0, 3.0, -3.0, 3.0, -3.0, 
3.0, -3.0), and X~ = (-1.5, 1.5,-1.5, 1.5, -1.5, 1.5, 
-1.5, 1.5, -1.5, 1.5) are stable and X° = (0, 0, 0, 0, 0, 



0, 0, 0, 0, 0) is unstable. For example considering Xp~, 
the Jacobian I|lj4(l takes the form 



J{Xj 



(- 


1.5 - D 


D 











1 











o\ 




D 


-1.5 - 2D 


D 











1 
















D 


-1.5 - 2D 


D 











1 
















D 


-1.5 -2D 


D 











1 
















D 


-1.5 - D 














1 




-1 














-1 



















-1 














-1 



















-1 














-1 



















-1 














-1 





\ 














-1 














-lj 



TABLE II: Stable steady states Xf and X~ for < D < 0.4703. Only x is ,i = 1,2... 5 are given. The values of y's are 
obtained from the relation y is — —x is (One may note that the second initial conditions {2:1(0), £2(0), 2:3(0), £4(0), 2:5(0)} 
= {2.7990, 2.1709, 1.7803, -1.4433, -1.4923} and yi(0) = -Xi(Q), i = 1,2, ..5, chosen in Sec. IIII B 1HI B lbl is closer to Xf 
than X~ . Similarly the first of the initial conditions (2:1(0), 2:2(0), 2:3(0), 2:4(0), 2:5(0)} = {3.0, 3.0, —1.5, —1.5, —1.5} and 
J/i(0) = — Xi(0), i = 1, 2. . . 5, is nearer to X~ than Xf.) 



D 




X2s 


2:3s 


.X'4 


S 


x$ 


3 


D 


Xls 


2:2s 


x 3 


3 


X4 


S 


Xs 


S 


Xf 0.0000 


3.0000 


3.0000 


3.0000 


-1 


5000 


-1 


5000 


X~ 0.0000 


3.0000 


3.0000 


-1 


5000 


-1 


5000 


-1 


5000 


0.0100 


3.0000 


2.9997 


2.9647 


-1 


4823 


-1 


4999 


0.0100 


2.9997 


2.9647 


-1 


4823 


-1 


4999 


-1 


5000 


0.0200 


3.0000 


2.9989 


2.9308 


-1 


4651 


-1 


4997 


0.0200 


2.9989 


2.9308 


-1 


4651 


-1 


4997 


-1 


5000 


0.0300 


2.9999 


2.9977 


2.8981 


-1 


4485 


-1 


4994 


0.0300 


2.9976 


2.8981 


-1 


4485 


-1 


4994 


-1 


5000 


0.0400 


2.9999 


2.9960 


2.8666 


-1 


4323 


-1 


4989 


0.0400 


2.9959 


2.8666 


-1 


4323 


-1 


4989 


-1 


5000 


0.0500 


2.9998 


2.9939 


2.8362 


-1 


4166 


-1 


4984 


0.0500 


2.9937 


2.8362 


-1 


4166 


-1 


4984 


-1 


5000 


0.4100 


2.9571 


2.8262 


2.1655 


-1 


0393 


-1 


4351 


0.4100 


2.7921 


2.1584 


-1 


0411 


-1 


4423 


-1 


4919 


0.4200 


2.9550 


2.8209 


2.1537 


-1 


0321 


-1 


4327 


0.4200 


2.7853 


2.1462 


-1 


0340 


-1 


4403 


-1 


4914 


0.4300 


2.9528 


2.8156 


2.1422 


-1 


0250 


-1 


4303 


0.4300 


2.7784 


2.1342 


-1 


0270 


-1 


4383 


-1 


4909 


0.4400 


2.9506 


2.8102 


2.1308 


-1 


0180 


-1 


4279 


0.4400 


2.7715 


2.1223 


-1 


0202 


-1 


4363 


-1 


4905 


0.4500 


2.9484 


2.8049 


2.1196 


-1 


0111 


-1 


4254 


0.4500 


2.7646 


2.1107 


-1 


0134 


-1 


4343 


-1 


4900 


0.4600 


2.9461 


2.7996 


2.1087 


-1 


0043 


-1 


4230 


0.4600 


2.7577 


2.0993 


-1 


0068 


-1 


4322 


-1 


4895 


0.4700 


2.9435 


2.7934 


2.0937 


-1 


0162 


-1 


4234 


0.4700 


2.7508 


2.0880 


-1 


0002 


-1 


4302 


-1 


4890 



The eigenvalues of the above matrix l|B5|) as a func- 
tion of D can be obtained by solving the corresponding 
characteristic equation. We have verified by numerical 
evaluation that all the eigenvalues of i|B5|) have negative 
real parts for D > 0. Similarly all the eigenvalues corre- 
sponding to the steady state X~ have negative real parts 
irrespective of the values of D (D > 0). We have also ver- 
ified that atleast one of the eigenvalues corresponding to 
the state X® has positive real part for D > 0. Thus Xp 
are stable, while X® is unstable for all values of D > 0. 



Similarly one can obtain the stability properties of the 
set of steady states X s — {Xf, X", X~}, by substituting 
the appropriate expressions from |Bj) into l|B4|) . Numer- 
ical diagonalization of i|B4|) for these steady states show 
that X® is unstable for all values of D, while Xf and 
X~ are stable for < D < 0.4703 and unstable for 
D > 0.4703 for the chosen initial conditions. The nu- 
merical values of Xjf- in the stable region are given in 
Table ITT1 as a function of D. 

x 
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